rm(list=ls())
library(tidyverse)
setwd("~/Downloads/dataverse_files(3)")

getmode <- function(v, na.omit = TRUE){
  if (na.omit == TRUE) {
    v <- na.omit(v)  
  }
  uv <- unique(v)
  tab <- tabulate(match(v, uv))
  out <- uv[tab == max(tab)]
  if(length(out) > 1){
    out <- sample(out,1)
  }
  return(out)
}

# add custom functions to extract estimates (tidy) and goodness-of-fit (glance) information
tidy.feglm <- function(x, ...) {
  s <- summary(x,  type = "sandwich")
  ret <- data.frame(
    term      = row.names(s$cm),
    estimate  = s$cm[, 1],
    std.error  = s$cm[, 2]
  )
  ret
}

glance.feglm <- function(x, ...) {
  ret <- data.frame(
    Null.Dev = x$null.deviance,
    Deviance = x$deviance,
    N   = x$nobs[4],
    `Committee FE` = as.integer(length(x$nms.fe[[1]])),
    `Congress FE` = as.integer(length(x$nms.fe[[2]])))
  ret
}


documents_supp <- read_csv("Data/CongressDocumentsAnalysisSupplementaryData.csv")
hearings<- read_csv("Data/CongressDocumentsCitations.csv")


hearings <- hearings %>% mutate(cong_doc_type = case_when(grepl("event", hearings$policy_document_url) == T ~ "Hearing",
                                                          grepl("report", hearings$policy_document_url) == T ~ "Report",
                                                          grepl("print", hearings$policy_document_url) == T ~ "Publication", 
                                                          TRUE ~ NA))

hearings %>% filter(cong_doc_type == "Hearing", congress > 106, congress < 114) %>% pull(policy_document_id) %>% unique() %>% length()

documents_supp %>% pull(policy_document_id) %>% unique() %>% length()
documents_supp$party_control <- fct_relevel(documents_supp$party_control, "R", "D")



library(fixest)
library(marginaleffects)


documents_supp$chair_extremism <- abs(documents_supp$median_floor - documents_supp$chair_dw)


doc_cites_model_minimal <- feglm(cites_science ~ party_control | committee_short + congress , data = documents_supp, family = binomial(link = "logit"), se = "twoway")

doc_cites_model_subcommittee <- feglm(cites_science ~ party_control | committee_short + congress, data = filter(documents_supp, subcommittee == 1), family = binomial(link = "logit"),  se = "twoway")

doc_cites_model_fullcommittee <- feglm(cites_science ~ party_control | committee_short  + congress, data = filter(documents_supp, subcommittee == 0), family = binomial(link = "logit"),  se = "twoway")

doc_cites_model_subcommittee_interaction <- feglm(cites_science ~ party_control*subcommittee | committee_short +congress, data = documents_supp, family = binomial(link = "logit"), se = "twoway")

library(modelsummary)
modelsummary(list("All Matched" = doc_cites_model_minimal, "Subcommittee" = doc_cites_model_subcommittee, "Full Committee" = doc_cites_model_fullcommittee, "Subommittee Interaction" = doc_cites_model_subcommittee_interaction), output = "Output/TableS22.docx", stars = T)

politicization <- read_csv("Data/ZetaBasedPoliticizationScores.csv")
voteunity <- read_csv("Data/VoteUnityScores.csv")


voteunity <- voteunity %>% group_by(pap_majortopic, cong) %>% summarise(topic_vote_unity = mean(party_unity, na.rm = T))

documents_supp <- documents_supp %>% left_join(dplyr::select(politicization, Major, congress, mean_wtd_abs_zeta) , by = c("majortopic" = "Major", "congress" = "congress"))

documents_supp <- documents_supp %>% left_join(voteunity, by = c("majortopic" = "pap_majortopic", "congress" = "cong"))


doc_cites_model_minimal <- feglm(cites_science ~ party_control | committee_short + congress , data = documents_supp, family = binomial(link = "logit"), se = "twoway")

doc_cites_model_vote_unity<- feglm(cites_science ~ party_control + topic_vote_unity | committee_short + congress  , data = documents_supp, family = binomial(link = "logit"), se = "twoway")
doc_cites_model_vote_unity_topicFE <- feglm(cites_science ~ party_control + topic_vote_unity | committee_short +  majortopic + congress , data = documents_supp, family = binomial(link = "logit"), se = "threeway")
doc_cites_model_vote_unity_int_topicFE <- feglm(cites_science ~ party_control*topic_vote_unity | committee_short + majortopic + congress , data = documents_supp, family = binomial(link = "logit"), se = "threeway")

doc_cites_model_textshare <- feglm(cites_science ~ party_control + mean_wtd_abs_zeta |committee_short +  congress  , data = documents_supp, family = binomial(link = "logit"), se = "twoway")
doc_cites_model_textshare_topicFE <- feglm(cites_science ~ party_control + mean_wtd_abs_zeta |committee_short + majortopic + congress  , data = documents_supp, family = binomial(link = "logit"), se = "threeway")
doc_cites_model_textshare_int_topicFE <- feglm(cites_science ~ party_control*mean_wtd_abs_zeta |committee_short + majortopic + congress  , data = documents_supp, family = binomial(link = "logit"), se = "threeway")



library(modelsummary)
modelsummary(models = list(doc_cites_model_minimal, doc_cites_model_vote_unity, doc_cites_model_vote_unity_topicFE, doc_cites_model_vote_unity_int_topicFE, 
                           doc_cites_model_textshare,doc_cites_model_textshare_topicFE, doc_cites_model_textshare_int_topicFE), output = "Output/TableS32.docx", stars = T)



zeta_pred <- feols(mean_wtd_abs_zeta ~ party_control | committee_short + congress + majortopic, data = documents_supp, se = "threeway")

summary(zeta_pred)
modelsummary(models = list(zeta_pred), output = "Output/TableS33.docx", stars = T)


doc_cites_model_allwitnessesD <- feglm(cites_science ~  n_agriculture + n_corporation + n_tradeassoc +      
                                         n_bureaucrat + n_congressional + n_judicial + n_statelocgov + n_education +        
                                         n_research + n_university + n_lawyerlobbyist + n_labor +  n_nonprofit +        
                                         n_healthcare + n_membership + n_citizen +  n_native +   n_religious +         
                                         n_international |congress +committee_short, data = filter(documents_supp, party_control == "D"), family = binomial(link = "logit"))


doc_cites_model_allwitnessesR <- feglm(cites_science ~    n_agriculture + n_corporation + n_tradeassoc +      
                                         n_bureaucrat + n_congressional + n_judicial + n_statelocgov + n_education +        
                                         n_research + n_university + n_lawyerlobbyist + n_labor +  n_nonprofit +        
                                         n_healthcare + n_membership + n_citizen +  n_native +   n_religious +         
                                         n_international | congress +committee_short , data = filter(documents_supp, party_control == "R"), family = binomial(link = "logit"))


witness_models <- modelplot(list("Democratic" = doc_cites_model_allwitnessesD, "Republican" = doc_cites_model_allwitnessesR), coef_omit = "(Intercept)", exponentiate = T) + 
  geom_vline(xintercept = 1, linetype = 3) + scale_color_manual("", values=c("#156B90","#9A3E25")) + theme(legend.position = "top") + xlab("Odds Ratio \n(95% CIs)")

ggsave("Output/FigS29.pdf", witness_models, width = 5, height = 6)



############ Partisan differences in sources in hearings
library(readstata13)

ban_dat <- readstata13::read.dta13("Data/hearings_apsr.dta")
ban_dat <- ban_dat %>% as_tibble()

ban_dat <- ban_dat %>% as_tibble() %>% filter(congress <117, congress > 106)
ban_dat <- ban_dat %>% filter(title != "")

ban_dat <- ban_dat %>% filter(congress > 106, congress < 115)
library(fixest)

n_university_mod <-  fenegbin(n_university ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_agriculture_mod <-  fenegbin(n_agriculture ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_tradeassoc_mod <-  fenegbin(n_tradeassoc ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_bureaucrat_mod <-  fenegbin(n_bureaucrat ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_congressional_mod <-  fenegbin(n_congressional ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_judicial_mod <-  fenegbin(n_judicial ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_statelocgov_mod <-  fenegbin(n_statelocgov ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_education_mod <-  fenegbin(n_education ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_research_mod <-  fenegbin(n_research ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_lawyerlobbyist_mod <-  fenegbin(n_lawyerlobbyist ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_labor_mod <-  fenegbin(n_labor ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_nonprofit_mod <-  fenegbin(n_nonprofit ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_healthcare_mod <-  fenegbin(n_healthcare ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_membership_mod <-  fenegbin(n_membership ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_citizen_mod <-  fenegbin(n_citizen ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_native_mod <-  fenegbin(n_native ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_religious_mod <-  fenegbin(n_religious ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")
n_international_mod <-  fenegbin(n_international ~ dem_maj | congress + cmt_thomas_id + majortopic, data = ban_dat, se = "threeway")



matching_objects <- grep("^n_.*_mod$", ls(), value = TRUE)
matching_objects <- mget(matching_objects, envir = .GlobalEnv)
tidy_models <- matching_objects %>% map(~tidy(.))
tidy_models <- tidy_models %>% bind_rows(.id = "model")

tidy_models$model <- gsub("_mod", "", tidy_models$model)
tidy_models$model <- gsub("n_", "", tidy_models$model)

tidy_models <- tidy_models %>% filter(term == "dem_maj") %>% mutate(model = fct_reorder(toupper(model), estimate))
tidy_models$sig <- 0
tidy_models$sig[tidy_models$p.value <= .05] <- 1

witness_type <- ggplot(tidy_models, aes(x=model, y = estimate, ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error, color = as.factor(sig))) +
  geom_pointrange() + coord_flip() + theme_minimal() + 
  geom_hline(yintercept = 0, linetype = 3) + scale_color_manual(values = c("darkgrey", "black")) +
  theme(legend.position = "none") + ylab("Coefficient Estimate for Democratic Control of Committee") + xlab("Witness Type")

ggsave("Output/FigS28.pdf", witness_type, width = 6, height = 5)

